%%

X1 = X_ref + 2*rand(length(X_ref),1)/115600000;
X1 = X1 / norm(X1);

%[X1, L2, rho] = powermethodRHO(A, X1, 9);
rho = 0.99403466159856/0.999999922699637;
%B   = -rho;
B   = -0.53105945;

%%

l = 1;

[a, b] = alpha_beta(rho, B, 1); 

V = A*X1; 
L2 = (V'*X1)/(X1'*X1); 
X2 = X1 + a*(V/L2 - X1) + b*(X1 - X0);

R1 = V/L2 - X1; 
norm(R1)

%%

l = 2;

X0 = X1; X1 = X2; L1 = L2; R1 = R2;

V  = A*X1; 
L2 = (V'*X1)/(X1'*X1); 
X2 = X1 + a*(V/L2 - X1) + b*(X1 - X0);

R2 = V/L2 - X1; 
norm(R2)

q = norm(R2)/norm(R1)
r = rho_est(rho, q, l)

%%

l = 3;

X0 = X1; X1 = X2; L1 = L2; R1 = R2;

V  = A*X1; 
L2 = (V'*X1)/(X1'*X1); 
X2 = X1 + a*(V/L2 - X1) + b*(X1 - X0);

R2 = V/L2 - X1; 
norm(R2)

q = norm(R2)/norm(R1)
r = rho_est(rho, q, l)

%%

l = 4;

X0 = X1; X1 = X2; L1 = L2; R1 = R2;

V  = A*X1; 
L2 = (V'*X1)/(X1'*X1); 
X2 = X1 + a*(V/L2 - X1) + b*(X1 - X0);

R2 = V/L2 - X1; 
norm(R2)

q = norm(R2)/norm(R1)
r = rho_est(rho, q, l)

%%

l = 5;

X0 = X1; X1 = X2; L1 = L2; R1 = R2;

V  = A*X1; 
L2 = (V'*X1)/(X1'*X1); 
X2 = X1 + a*(V/L2 - X1) + b*(X1 - X0);

R2 = V/L2 - X1; 
norm(R2)

q = norm(R2)/norm(R1)
r = rho_est(rho, q, l)

%%

l = 6;

X0 = X1; X1 = X2; L1 = L2; R1 = R2;

V  = A*X1; 
L2 = (V'*X1)/(X1'*X1); 
X2 = X1 + a*(V/L2 - X1) + b*(X1 - X0);

R2 = V/L2 - X1; 
norm(R2)

q = norm(R2)/norm(R1)
r = rho_est(rho, q, l)